R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/define_heatwaves.m

Note: “00_OISST_include_plotx.R” contains included library and gplotx function in previous Rmd files.

yrng <- seq(1982,2022)
clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))

#### Define Marine heatwave is monthly anomaly > 90% historical anomaly
Thresh <- 0.9
Detrend <- TRUE #FALSE
if (Detrend) {
  prex <- "../data_src/oisst/monthly_anom_detrend/"
  outd <- "../data_src/oisst/monthly_heatwave_detrend/"
} else {
  prex <- "../data_src/oisst/monthly_anom_icemask/"
  outd <- "../data_src/oisst/monthly_heatwave/"
}

Recalculate_thre <- FALSE
if (Recalculate_thre) {
  
# It's a 3-month running window, and use quantile(probs=Thresh=0.9) to aggregate the window and get the threshold(percentile>0.9) value
sty <- future_lapply(1:12, function(j) {
  if (j==1) {
    winj = c(12, 1, 2)
  } else if (j==12) {
    winj = c(11, 12, 1)
  } else {
    winj = c(j-1, j , j+1)
  }
  monj <- fifelse(winj<10, paste0("0",winj), paste0(winj))
  jstr <- monstr[j]
  for (i in clim_years) {
    if (winj[1] == 12 & (i-1)>=clim_years[1]) {
      x0 <- read_stars(paste0(prex, i-1, monj[1], "_anom.nc"))
    } else if (winj[1] == 12 & i==clim_years[1]) {
      x0 <- NULL
    } else {
      x0 <- read_stars(paste0(prex, i, monj[1], "_anom.nc"))
    }
    
    x1 <- read_stars(paste0(prex, i, monj[2], "_anom.nc"))
    
    if (winj[3] == 1 & (i+1)<=clim_years[length(clim_years)]) {
      x2 <- read_stars(paste0(prex, i+1, monj[3], "_anom.nc"))
    } else if (winj[3] == 1 & i==clim_years[length(clim_years)]) {
      x2 <- NULL
    } else {
      x2 <- read_stars(paste0(prex, i, monj[3], "_anom.nc"))
    }
    
    if (is.null(x0)) {
      x <- c(x1, x2)
      datex<- c(as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"),
                as.Date(paste0(i, monj[3], "01"), format="%Y%m%d"))
    } else if (is.null(x2)) {
      x <- c(x0, x1)
      datex<- c(as.Date(paste0(i, monj[1], "01"), format="%Y%m%d"),
                as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"))
    } else {
      x <- c(x0, x1, x2)
      datex<- c(fifelse(winj[1]==12, as.Date(paste0(i-1, monj[1], "01"), format="%Y%m%d"), as.Date(paste0(i, monj[1], "01"), format="%Y%m%d")),
                as.Date(paste0(i, monj[2], "01"), format="%Y%m%d"),
                fifelse(winj[3]==1, as.Date(paste0(i+1, monj[3], "01"), format="%Y%m%d"), as.Date(paste0(i, monj[3], "01"), format="%Y%m%d")))
    }
    if (i == clim_years[1]) {
      styx <- x
      datey<- datex
    } else {
      styx <- c(styx, x)
      datey<- c(datey, datex)
    }
  }
  print(paste0("Before process quantile, month: ", monj, " ..."))
  
  names(styx) <- rep(jstr, length(names(styx)))
  styx <- merge(styx) %>% st_set_dimensions(3, values = as.POSIXct(datey), names = "time") %>% 
    aggregate(by=paste0(climyrs, " years"), FUN=quantile, na.rm=TRUE, probs=Thresh, names=FALSE)
  names(styx)[1] <- jstr
  styx <- styx %>% dplyr::select(jstr) %>% adrop
  return (styx)
})

  save(sty, file="../data_src/stats/sst_1982_2011mhw_threshold_nc.RData")
} else {
  load("../data_src/stats/sst_1982_2011mhw_threshold_nc.RData")
}
Recalculate_mhw <- FALSE

if (Recalculate_mhw) {
for (i in yrng) {
  mmx <- fifelse(i==curryr, currmo-1L, 12L)
  for (j in 1:mmx) {
    monj <- fifelse(j<10, paste0("0",j), paste0(j))
    jstr <- monstr[j]
    print(paste0("Now in Year-month: ", i," - ", monj, " to get heatwave"))
    z <- read_stars(paste0(prex, i, monj, "_anom.nc"))
    names(z)[1] <- "heatwave" 
    zt <- matrix(rep(0, len=1440*720), nrow = 1440)
    zt[which(z[[1]]>=sty[[j]][[1]] & z[[1]]>0)] <- 1 ## NOTE: HERE different with ref code, add criteria: ANOMALY > 0 (for MHWs)
    z[[1]] <- zt

    filet <- paste0(outd, i, monj, "_heatwave.nc")
    write_stars(z, filet)
  }
}
}
ht <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/202204_heatwave.nc"))
names(ht)[1] <- "heatwave"

dht <- as.data.table(ht) %>% 
  .[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, heatwave)]

ghx <- gplotx(dht, "heatwave", returnx=TRUE, minz=-0.4, maxz=1, maxlim=180, no_coord=TRUE, legend_label="2022-04 Marine Heatwaves  ", legend_pos=c(0.65, 0.09))
ghx <- ghx + geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3)

ghx

zk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "02", "_anom.nc"))
zt <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "02", "_anom.nc"))
names(zk)[1] <- "anom"
names(zt)[1] <- "anom"

zkt <- as.data.table(zk) %>% .[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, anom)]
ztt <- as.data.table(zt) %>% .[,`:=`(x=fifelse(x>180, x-360, x))] %>% .[,.(x, y, anom)]

#print(range(na.omit(zkt$anom))) #-9.520356  7.962329
#print(range(na.omit(ztt$anom))) #-5.414660  7.204932

gt <- gplotx(ztt, "anom", returnx=TRUE, minz=-5.5, maxz=7.5, fillopts=NA, legend_label="2020-02 Anomaly  ", legend_pos=c(0.65, 0.09), maxlim=180)

gk <- gplotx(zkt, "anom", returnx=TRUE, minz=-10, maxz=7.5, fillopts=NA, legend_label="2020-02 Detrended Anomaly  ", legend_pos=c(0.65, 0.09), maxlim=180)

layt <- rbind(c(1,1),c(2,2),c(3,3))

grid.arrange(gt, gk, ghx,  layout_matrix=layt)

if (!require("layer")) {
  if (!require("remotes")) install.packages("remotes"); library(remotes)
  remotes::install_github("marcosci/layer")
  library(layer)
}

Reload_hw_nc <- FALSE
if (Reload_hw_nc) {
  hw202002 <-  as.data.table(ht) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    #.[,.(longitude, latitude, heatwave)] %>%
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  ht1 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/201002_heatwave.nc"))
  names(ht1)[1] <- "heatwave"

  hw201002 <-  as.data.table(ht1) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  ht2 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/200002_heatwave.nc"))
  names(ht2)[1] <- "heatwave"

  hw200002 <-  as.data.table(ht2) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  ht3 <- read_stars(paste0("../data_src/oisst/monthly_heatwave_detrend/199002_heatwave.nc"))
  names(ht3)[1] <- "heatwave"

  hw199002 <-  as.data.table(ht3) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))
  
  save(hw199002, hw200002, hw201002, hw202002, file="../data_src/stats/heatwave_1990_2020_stack.RData")
} else {
  load("../data_src/stats/heatwave_1990_2020_stack.RData")
}

tl0 <- layer::tilt_map(hw199002)
tl1 <- layer::tilt_map(hw200002, y_shift=75)
tl2 <- layer::tilt_map(hw201002, y_shift=150)
tl3 <- layer::tilt_map(hw202002, y_shift=225)

map_list <- list(tl0, tl1, tl2, tl3)
#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "cividis")